library(SingleCellExperiment)
library(scater)
library(scran)
library(BiocParallel)
sce <- readRDS("data/sce_annotated.rds")Subset to T cells
sce_big <- sce
sce <- sce[, sce$label_immgen == "T cells"]
sce## class: SingleCellExperiment
## dim: 32286 3225
## metadata(1): Samples
## assays(2): counts logcounts
## rownames(32286): ENSMUSG00000051951 ENSMUSG00000089699 ...
## ENSMUSG00000095041 YFP
## rowData names(3): ID Symbol Type
## colnames: NULL
## colData names(6): Sample Barcode ... label_immgen
## label_immgen_collapsed
## reducedDimNames(2): PCA UMAP
## mainExpName: NULL
## altExpNames(0):Rerun UMAP and clustering
#--- variance-modelling ---#
set.seed(00010101)
dec.sce <- modelGeneVarByPoisson(sce)
top.sce <- getTopHVGs(dec.sce, prop=0.1)
set.seed(101010011)
sce <- denoisePCA(sce, technical=dec.sce, subset.row=top.sce)
#sce <- runTSNE(sce, dimred="PCA")
sce <- runUMAP(sce, dimred="PCA",
metric = "cosine",
min_dist = 0.3,
n_neighbors = 30)
snn.gr <- buildSNNGraph(sce, use.dimred="PCA", k=25)
set.seed(2)
colLabels(sce) <- factor(igraph::cluster_leiden(snn.gr, resolution_parameter = 1.8)$membership)
table(colLabels(sce))##
## 1 2 3 4 5 6 7 8 9
## 542 795 521 191 469 88 510 105 4sce <- sce[, sce$label != 9]
sce$label <- droplevels(sce$label)Define custom colors:
cluster_colors = c(
`1` = '#fdbf6f',
`2` = '#b2df8a',
`3` = '#33a02c',
`4` = '#fb9a99',
`5` = '#e31a1c',
`6` = '#a6cee3',
`7` = '#1f78b4',
`8` = '#ff7f00'
)plotUMAP(sce, colour_by="label", text_by="label") + coord_fixed() +
ggh4x::force_panelsizes(rows = unit(4, "in"),
cols = unit(4, "in")) +
scale_color_manual(values = cluster_colors)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plotUMAP(sce, colour_by="label", other_fields="Sample") + coord_fixed() +
facet_wrap(~Sample) +
ggh4x::force_panelsizes(rows = unit(2, "in"),
cols = unit(2, "in")) +
scale_color_manual(values = cluster_colors)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plotTableBarplot <- function(x, y, label_x, label_y) {
tbl <- table(x, y)
m <- prop.table(tbl, margin = 1)*100
m <- as.data.frame(m)
ggplot(m, aes(x=x, y=Freq, fill=y)) +
geom_bar(stat = "identity",
color="#444444") +
xlab(label_x) +
ylab(paste( "%" , label_y)) +
guides(fill=guide_legend(title=label_y)) +
ggthemes::theme_few() +
#scale_fill_manual(values = scater:::.get_palette("tableau10medium"))
scale_fill_manual(values = cluster_colors)
}
plotTableBarplot(x = sce$Sample,
y = sce$label,
label_x = "Sample",
label_y = "Cluster") +
ggh4x::force_panelsizes(rows = unit(4, "in"),
cols = unit(2, "in")) +
xlab('')
MAGIC
library(magieR)
sce <- magieR(sce, n.jobs=14)We copy the rowData, colData and UMAP to the magic dataset
rowData(altExp(sce, "magic")) <- rowData(sce)
colData(altExp(sce, "magic")) <- colData(sce)
reducedDim(altExp(sce, "magic"), "UMAP") <- reducedDim(sce, "UMAP")Marker genes
Plot some marker genes to get an idea about the cell types
library(patchwork)
genes <- c("Cd4" , "Cd8a", "Foxp3", "Pdcd1", "Tcf7", "Adrb1", "Adrb2",
"Havcr2", "Cd101", "Cx3cr1", "Cxcr6", "Slamf6", "Entpd1", "Gata3", "Il2ra", "Mki67",
"Tox", "Cxcr3", "Gzmb", "Tbx21", "Cd44", "Klrk1", "Prf1", "Lag3",
"Sell", "Itgae", "Klrg1", "Bcl6", "Cxcr5", "Tnfrsf4", "Cd200", "Itga1",
"Ifng", "Tnf", "Gzma", "Gzmb", "Prf1", "Il2", 'Ikzf2', "Klf3", "Runx3",
"Klf2", "Nr4a1", "S1pr1", "Cd69", "Zfp683", "Prdm1")
genes = unique(genes)
genes = sort(genes)
plots <- purrr::map(genes, ~ {
p_umap <- plotUMAP(altExp(sce),
swap_rownames = "Symbol",
colour_by = .x) +
ggtitle(.x) +
ggthemes::theme_few() +
coord_fixed() +
theme(legend.position = c(0.8, 0.2),
#legend.spacing.x = unit(0, "mm"),
legend.spacing.y = unit(0, "mm"),
legend.title = element_blank(),
legend.background = element_rect(fill=alpha('white', 0.4)),
legend.text = element_text(size=8),
plot.title = element_text(face = "bold.italic"))
#p_umap
p_violin <- plotExpression(
altExp(sce),
features = .x,
swap_rownames = "Symbol",
colour_by = "label",
x = "label"
) +
ggthemes::theme_few() +
xlab("Cluster") +
ylab("Expression") +
theme(strip.background = element_blank(),
strip.text = element_blank(),
legend.position = "none") +
scale_color_manual(values = cluster_colors)
list(p_umap, p_violin)
})
plots <- unlist(plots, recursive = FALSE)
plot_cluster <- plotUMAP(sce, colour_by="label", text_by="label") + ggthemes::theme_few() +
ggtitle("Cluster") +
theme(legend.position = "none",
plot.title = element_text(face = "bold")) +
scale_color_manual(values = cluster_colors) +
coord_fixed()
plots[[length(plots) + 1]] = plot_cluster
p <- patchwork::wrap_plots(plots, ncol = 8) & theme(legend.key.size = unit(2.5, 'mm'))
p
Heatmap of DE-genes
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:Biobase':
##
## combine## The following objects are masked from 'package:GenomicRanges':
##
## intersect, setdiff, union## The following object is masked from 'package:GenomeInfoDb':
##
## intersect## The following objects are masked from 'package:IRanges':
##
## collapse, desc, intersect, setdiff, slice, union## The following objects are masked from 'package:S4Vectors':
##
## first, intersect, rename, setdiff, setequal, union## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, union## The following object is masked from 'package:matrixStats':
##
## count## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ComplexHeatmap)
markers.up <- purrr::map(unique(colLabels(sce)), ~{
res = findMarkers(sce,
groups = colLabels(sce) == .x,
pval.type="all", direction="up", row.data=rowData(sce),
lfc = 0.1,
add.summary = TRUE,
test.type = "binom")
res$`TRUE`
})
markers.up## [[1]]
## DataFrame with 32286 rows and 11 columns
## ID Symbol Type self.average
## <character> <character> <character> <numeric>
## ENSMUSG00000035042 ENSMUSG00000035042 Ccl5 Gene Expression 4.950571
## ENSMUSG00000004612 ENSMUSG00000004612 Nkg7 Gene Expression 2.695000
## ENSMUSG00000042385 ENSMUSG00000042385 Gzmk Gene Expression 0.343448
## ENSMUSG00000053977 ENSMUSG00000053977 Cd8a Gene Expression 1.722286
## ENSMUSG00000048521 ENSMUSG00000048521 Cxcr6 Gene Expression 0.885217
## ... ... ... ... ...
## ENSMUSG00000095763 ENSMUSG00000095763 AC124606.2 Gene Expression 0
## ENSMUSG00000095523 ENSMUSG00000095523 AC124606.1 Gene Expression 0
## ENSMUSG00000095475 ENSMUSG00000095475 AC133095.2 Gene Expression 0
## ENSMUSG00000094855 ENSMUSG00000094855 AC133095.1 Gene Expression 0
## ENSMUSG00000095019 ENSMUSG00000095019 AC234645.1 Gene Expression 0
## other.average self.detected other.detected p.value
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000035042 0.9524074 0.970480 0.2952594 7.48249e-79
## ENSMUSG00000004612 0.9748792 0.837638 0.3613289 6.59609e-38
## ENSMUSG00000042385 0.0390605 0.145756 0.0238895 9.53273e-24
## ENSMUSG00000053977 0.7068299 0.649446 0.3184024 2.77388e-22
## ENSMUSG00000048521 0.2640694 0.357934 0.1414707 2.74553e-20
## ... ... ... ... ...
## ENSMUSG00000095763 0 0 0 1
## ENSMUSG00000095523 0 0 0 1
## ENSMUSG00000095475 0 0 0 1
## ENSMUSG00000094855 0 0 0 1
## ENSMUSG00000095019 0 0 0 1
## FDR summary.logFC logFC.FALSE
## <numeric> <numeric> <numeric>
## ENSMUSG00000035042 2.41580e-74 1.71461 1.71461
## ENSMUSG00000004612 1.06481e-33 1.21161 1.21161
## ENSMUSG00000042385 1.02591e-19 2.57822 2.57822
## ENSMUSG00000053977 2.23894e-18 1.02693 1.02693
## ENSMUSG00000048521 1.77285e-16 1.33537 1.33537
## ... ... ... ...
## ENSMUSG00000095763 1 0 0
## ENSMUSG00000095523 1 0 0
## ENSMUSG00000095475 1 0 0
## ENSMUSG00000094855 1 0 0
## ENSMUSG00000095019 1 0 0
##
## [[2]]
## DataFrame with 32286 rows and 11 columns
## ID Symbol Type self.average
## <character> <character> <character> <numeric>
## ENSMUSG00000023927 ENSMUSG00000023927 Satb1 Gene Expression 1.937062
## ENSMUSG00000045092 ENSMUSG00000045092 S1pr1 Gene Expression 1.177296
## ENSMUSG00000026581 ENSMUSG00000026581 Sell Gene Expression 1.116014
## ENSMUSG00000027985 ENSMUSG00000027985 Lef1 Gene Expression 1.019217
## ENSMUSG00000020846 ENSMUSG00000020846 Rflnb Gene Expression 0.435921
## ... ... ... ... ...
## ENSMUSG00000095763 ENSMUSG00000095763 AC124606.2 Gene Expression 0
## ENSMUSG00000095523 ENSMUSG00000095523 AC124606.1 Gene Expression 0
## ENSMUSG00000095475 ENSMUSG00000095475 AC133095.2 Gene Expression 0
## ENSMUSG00000094855 ENSMUSG00000094855 AC133095.1 Gene Expression 0
## ENSMUSG00000095019 ENSMUSG00000095019 AC234645.1 Gene Expression 0
## other.average self.detected other.detected p.value
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000023927 0.735486 0.613836 0.3211047 1.23084e-22
## ENSMUSG00000045092 0.413271 0.411321 0.1908491 3.88949e-21
## ENSMUSG00000026581 0.481122 0.391195 0.2081616 2.24428e-14
## ENSMUSG00000027985 0.423690 0.364780 0.1986810 1.14798e-12
## ENSMUSG00000020846 0.128670 0.161006 0.0630668 1.30577e-12
## ... ... ... ... ...
## ENSMUSG00000095763 0 0 0 1
## ENSMUSG00000095523 0 0 0 1
## ENSMUSG00000095475 0 0 0 1
## ENSMUSG00000094855 0 0 0 1
## ENSMUSG00000095019 0 0 0 1
## FDR summary.logFC logFC.FALSE
## <numeric> <numeric> <numeric>
## ENSMUSG00000023927 3.97389e-18 0.933482 0.933482
## ENSMUSG00000045092 6.27880e-17 1.105321 1.105321
## ENSMUSG00000026581 2.41529e-10 0.908175 0.908175
## ENSMUSG00000027985 8.43159e-09 0.874524 0.874524
## ENSMUSG00000020846 8.43159e-09 1.343583 1.343583
## ... ... ... ...
## ENSMUSG00000095763 1 0 0
## ENSMUSG00000095523 1 0 0
## ENSMUSG00000095475 1 0 0
## ENSMUSG00000094855 1 0 0
## ENSMUSG00000095019 1 0 0
##
## [[3]]
## DataFrame with 32286 rows and 11 columns
## ID Symbol Type self.average
## <character> <character> <character> <numeric>
## ENSMUSG00000026581 ENSMUSG00000026581 Sell Gene Expression 1.54562
## ENSMUSG00000000782 ENSMUSG00000000782 Tcf7 Gene Expression 1.92293
## ENSMUSG00000027985 ENSMUSG00000027985 Lef1 Gene Expression 1.25667
## ENSMUSG00000023927 ENSMUSG00000023927 Satb1 Gene Expression 2.02189
## ENSMUSG00000026989 ENSMUSG00000026989 Dapl1 Gene Expression 0.72507
## ... ... ... ... ...
## ENSMUSG00000095763 ENSMUSG00000095763 AC124606.2 Gene Expression 0
## ENSMUSG00000095523 ENSMUSG00000095523 AC124606.1 Gene Expression 0
## ENSMUSG00000095475 ENSMUSG00000095475 AC133095.2 Gene Expression 0
## ENSMUSG00000094855 ENSMUSG00000094855 AC133095.1 Gene Expression 0
## ENSMUSG00000095019 ENSMUSG00000095019 AC234645.1 Gene Expression 0
## other.average self.detected other.detected p.value
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000026581 0.462654 0.608445 0.1848148 6.27945e-49
## ENSMUSG00000000782 0.724573 0.727447 0.2844444 4.49905e-39
## ENSMUSG00000027985 0.438306 0.518234 0.1859259 7.69630e-33
## ENSMUSG00000023927 0.841055 0.735125 0.3274074 1.92378e-30
## ENSMUSG00000026989 0.182960 0.295585 0.0766667 7.41670e-30
## ... ... ... ... ...
## ENSMUSG00000095763 0 0 0 1
## ENSMUSG00000095523 0 0 0 1
## ENSMUSG00000095475 0 0 0 1
## ENSMUSG00000094855 0 0 0 1
## ENSMUSG00000095019 0 0 0 1
## FDR summary.logFC logFC.FALSE
## <numeric> <numeric> <numeric>
## ENSMUSG00000026581 2.02738e-44 1.71568 1.71568
## ENSMUSG00000000782 7.26282e-35 1.35278 1.35278
## ENSMUSG00000027985 8.28275e-29 1.47579 1.47579
## ENSMUSG00000023927 1.55278e-26 1.16539 1.16539
## ENSMUSG00000026989 4.78911e-26 1.93829 1.93829
## ... ... ... ...
## ENSMUSG00000095763 1 0 0
## ENSMUSG00000095523 1 0 0
## ENSMUSG00000095475 1 0 0
## ENSMUSG00000094855 1 0 0
## ENSMUSG00000095019 1 0 0
##
## [[4]]
## DataFrame with 32286 rows and 11 columns
## ID Symbol Type
## <character> <character> <character>
## ENSMUSG00000022584 ENSMUSG00000022584 Ly6c2 Gene Expression
## ENSMUSG00000032666 ENSMUSG00000032666 1700025G04Rik Gene Expression
## ENSMUSG00000057329 ENSMUSG00000057329 Bcl2 Gene Expression
## ENSMUSG00000004612 ENSMUSG00000004612 Nkg7 Gene Expression
## ENSMUSG00000024910 ENSMUSG00000024910 Ctsw Gene Expression
## ... ... ... ...
## ENSMUSG00000095763 ENSMUSG00000095763 AC124606.2 Gene Expression
## ENSMUSG00000095523 ENSMUSG00000095523 AC124606.1 Gene Expression
## ENSMUSG00000095475 ENSMUSG00000095475 AC133095.2 Gene Expression
## ENSMUSG00000094855 ENSMUSG00000094855 AC133095.1 Gene Expression
## ENSMUSG00000095019 ENSMUSG00000095019 AC234645.1 Gene Expression
## self.average other.average self.detected other.detected
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000022584 1.552598 0.294772 0.596859 0.117822
## ENSMUSG00000032666 0.318471 0.054631 0.193717 0.029703
## ENSMUSG00000057329 1.398871 0.634577 0.654450 0.274257
## ENSMUSG00000004612 2.232064 1.203323 0.858639 0.415182
## ENSMUSG00000024910 0.995101 0.446038 0.518325 0.217492
## ... ... ... ... ...
## ENSMUSG00000095763 0 0 0 0
## ENSMUSG00000095523 0 0 0 0
## ENSMUSG00000095475 0 0 0 0
## ENSMUSG00000094855 0 0 0 0
## ENSMUSG00000095019 0 0 0 0
## p.value FDR summary.logFC logFC.FALSE
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000022584 1.41094e-35 4.55536e-31 2.33470 2.33470
## ENSMUSG00000032666 2.16899e-15 3.50141e-11 2.68004 2.68004
## ENSMUSG00000057329 2.63198e-14 2.83254e-10 1.25286 1.25286
## ENSMUSG00000004612 2.07575e-13 1.67544e-09 1.04720 1.04720
## ENSMUSG00000024910 1.20218e-11 7.76269e-08 1.25051 1.25051
## ... ... ... ... ...
## ENSMUSG00000095763 1 1 0 0
## ENSMUSG00000095523 1 1 0 0
## ENSMUSG00000095475 1 1 0 0
## ENSMUSG00000094855 1 1 0 0
## ENSMUSG00000095019 1 1 0 0
##
## [[5]]
## DataFrame with 32286 rows and 11 columns
## ID Symbol Type self.average
## <character> <character> <character> <numeric>
## ENSMUSG00000053310 ENSMUSG00000053310 Nrgn Gene Expression 1.231789
## ENSMUSG00000026475 ENSMUSG00000026475 Rgs16 Gene Expression 1.216841
## ENSMUSG00000030124 ENSMUSG00000030124 Lag3 Gene Expression 1.049661
## ENSMUSG00000026826 ENSMUSG00000026826 Nr4a2 Gene Expression 0.989286
## ENSMUSG00000041515 ENSMUSG00000041515 Irf8 Gene Expression 0.977445
## ... ... ... ... ...
## ENSMUSG00000095763 ENSMUSG00000095763 AC124606.2 Gene Expression 0
## ENSMUSG00000095523 ENSMUSG00000095523 AC124606.1 Gene Expression 0
## ENSMUSG00000095475 ENSMUSG00000095475 AC133095.2 Gene Expression 0
## ENSMUSG00000094855 ENSMUSG00000094855 AC133095.1 Gene Expression 0
## ENSMUSG00000095019 ENSMUSG00000095019 AC234645.1 Gene Expression 0
## other.average self.detected other.detected p.value
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000053310 0.064212 0.575693 0.0301599 1.80962e-143
## ENSMUSG00000026475 0.102185 0.590618 0.0483285 2.91712e-124
## ENSMUSG00000030124 0.132325 0.582090 0.0639535 8.45937e-106
## ENSMUSG00000026826 0.118990 0.509595 0.0559593 7.01253e-93
## ENSMUSG00000041515 0.137297 0.535181 0.0690407 1.26588e-88
## ... ... ... ... ...
## ENSMUSG00000095763 0 0 0 1
## ENSMUSG00000095523 0 0 0 1
## ENSMUSG00000095475 0 0 0 1
## ENSMUSG00000094855 0 0 0 1
## ENSMUSG00000095019 0 0 0 1
## FDR summary.logFC logFC.FALSE
## <numeric> <numeric> <numeric>
## ENSMUSG00000053310 5.84254e-139 4.22675 4.22675
## ENSMUSG00000026475 4.70910e-120 3.59438 3.59438
## ENSMUSG00000030124 9.10398e-102 3.17374 3.17374
## ENSMUSG00000026826 5.66017e-89 3.17274 3.17274
## ENSMUSG00000041515 8.17407e-85 2.94326 2.94326
## ... ... ... ...
## ENSMUSG00000095763 1 0 0
## ENSMUSG00000095523 1 0 0
## ENSMUSG00000095475 1 0 0
## ENSMUSG00000094855 1 0 0
## ENSMUSG00000095019 1 0 0
##
## [[6]]
## DataFrame with 32286 rows and 11 columns
## ID Symbol Type self.average
## <character> <character> <character> <numeric>
## ENSMUSG00000029075 ENSMUSG00000029075 Tnfrsf4 Gene Expression 1.869825
## ENSMUSG00000026770 ENSMUSG00000026770 Il2ra Gene Expression 0.977577
## ENSMUSG00000040204 ENSMUSG00000040204 Pclaf Gene Expression 0.889277
## ENSMUSG00000037706 ENSMUSG00000037706 Cd81 Gene Expression 0.783902
## ENSMUSG00000028832 ENSMUSG00000028832 Stmn1 Gene Expression 1.174229
## ... ... ... ... ...
## ENSMUSG00000095763 ENSMUSG00000095763 AC124606.2 Gene Expression 0
## ENSMUSG00000095523 ENSMUSG00000095523 AC124606.1 Gene Expression 0
## ENSMUSG00000095475 ENSMUSG00000095475 AC133095.2 Gene Expression 0
## ENSMUSG00000094855 ENSMUSG00000094855 AC133095.1 Gene Expression 0
## ENSMUSG00000095019 ENSMUSG00000095019 AC234645.1 Gene Expression 0
## other.average self.detected other.detected p.value
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000029075 0.316400 0.863636 0.1461858 1.46069e-29
## ENSMUSG00000026770 0.130271 0.568182 0.0683051 1.07582e-25
## ENSMUSG00000040204 0.137652 0.568182 0.0692627 1.84186e-25
## ENSMUSG00000037706 0.133870 0.522727 0.0683051 1.73813e-22
## ENSMUSG00000028832 0.248327 0.659091 0.1200128 1.35741e-21
## ... ... ... ... ...
## ENSMUSG00000095763 0 0 0 1
## ENSMUSG00000095523 0 0 0 1
## ENSMUSG00000095475 0 0 0 1
## ENSMUSG00000094855 0 0 0 1
## ENSMUSG00000095019 0 0 0 1
## FDR summary.logFC logFC.FALSE
## <numeric> <numeric> <numeric>
## ENSMUSG00000029075 4.71598e-25 2.55754 2.55754
## ENSMUSG00000026770 1.73669e-21 3.04481 3.04481
## ENSMUSG00000040204 1.98221e-21 3.02490 3.02490
## ENSMUSG00000037706 1.40293e-18 2.92465 2.92465
## ENSMUSG00000028832 8.76506e-18 2.45120 2.45120
## ... ... ... ...
## ENSMUSG00000095763 1 0 0
## ENSMUSG00000095523 1 0 0
## ENSMUSG00000095475 1 0 0
## ENSMUSG00000094855 1 0 0
## ENSMUSG00000095019 1 0 0
##
## [[7]]
## DataFrame with 32286 rows and 11 columns
## ID Symbol Type self.average
## <character> <character> <character> <numeric>
## ENSMUSG00000029075 ENSMUSG00000029075 Tnfrsf4 Gene Expression 1.197670
## ENSMUSG00000055435 ENSMUSG00000055435 Maf Gene Expression 1.013128
## ENSMUSG00000039521 ENSMUSG00000039521 Foxp3 Gene Expression 0.601526
## ENSMUSG00000025997 ENSMUSG00000025997 Ikzf2 Gene Expression 1.505087
## ENSMUSG00000002489 ENSMUSG00000002489 Tiam1 Gene Expression 0.340345
## ... ... ... ... ...
## ENSMUSG00000095763 ENSMUSG00000095763 AC124606.2 Gene Expression 0
## ENSMUSG00000095523 ENSMUSG00000095523 AC124606.1 Gene Expression 0
## ENSMUSG00000095475 ENSMUSG00000095475 AC133095.2 Gene Expression 0
## ENSMUSG00000094855 ENSMUSG00000094855 AC133095.1 Gene Expression 0
## ENSMUSG00000095019 ENSMUSG00000095019 AC234645.1 Gene Expression 0
## other.average self.detected other.detected p.value
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000029075 0.2010385 0.472549 0.1080782 1.28356e-52
## ENSMUSG00000055435 0.1822689 0.413725 0.0922169 1.70516e-47
## ENSMUSG00000039521 0.0532864 0.256863 0.0328292 4.19128e-46
## ENSMUSG00000025997 0.3271730 0.505882 0.1552932 7.91963e-40
## ENSMUSG00000002489 0.0430825 0.172549 0.0276651 3.28875e-27
## ... ... ... ... ...
## ENSMUSG00000095763 0 0 0 1
## ENSMUSG00000095523 0 0 0 1
## ENSMUSG00000095475 0 0 0 1
## ENSMUSG00000094855 0 0 0 1
## ENSMUSG00000095019 0 0 0 1
## FDR summary.logFC logFC.FALSE
## <numeric> <numeric> <numeric>
## ENSMUSG00000029075 4.14409e-48 2.12202 2.12202
## ENSMUSG00000055435 2.75264e-43 2.15805 2.15805
## ENSMUSG00000039521 4.51066e-42 2.94440 2.94440
## ENSMUSG00000025997 6.39233e-36 1.69982 1.69982
## ENSMUSG00000002489 2.12361e-23 2.61403 2.61403
## ... ... ... ...
## ENSMUSG00000095763 1 0 0
## ENSMUSG00000095523 1 0 0
## ENSMUSG00000095475 1 0 0
## ENSMUSG00000094855 1 0 0
## ENSMUSG00000095019 1 0 0
##
## [[8]]
## DataFrame with 32286 rows and 11 columns
## ID Symbol Type self.average
## <character> <character> <character> <numeric>
## ENSMUSG00000040204 ENSMUSG00000040204 Pclaf Gene Expression 1.76797
## ENSMUSG00000017716 ENSMUSG00000017716 Birc5 Gene Expression 1.19243
## ENSMUSG00000031004 ENSMUSG00000031004 Mki67 Gene Expression 1.63012
## ENSMUSG00000069272 ENSMUSG00000069272 Hist1h2ae Gene Expression 1.26500
## ENSMUSG00000094777 ENSMUSG00000094777 Hist1h2ap Gene Expression 1.77379
## ... ... ... ... ...
## ENSMUSG00000095763 ENSMUSG00000095763 AC124606.2 Gene Expression 0
## ENSMUSG00000095523 ENSMUSG00000095523 AC124606.1 Gene Expression 0
## ENSMUSG00000095475 ENSMUSG00000095475 AC133095.2 Gene Expression 0
## ENSMUSG00000094855 ENSMUSG00000094855 AC133095.1 Gene Expression 0
## ENSMUSG00000095019 ENSMUSG00000095019 AC234645.1 Gene Expression 0
## other.average self.detected other.detected p.value
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000040204 0.1039424 0.761905 0.0600128 1.15467e-50
## ENSMUSG00000017716 0.0685681 0.647619 0.0430039 4.50611e-47
## ENSMUSG00000031004 0.1325442 0.714286 0.0718870 1.36437e-41
## ENSMUSG00000069272 0.0537482 0.523810 0.0304878 1.06784e-40
## ENSMUSG00000094777 0.1016822 0.600000 0.0503851 5.98853e-39
## ... ... ... ... ...
## ENSMUSG00000095763 0 0 0 1
## ENSMUSG00000095523 0 0 0 1
## ENSMUSG00000095475 0 0 0 1
## ENSMUSG00000094855 0 0 0 1
## ENSMUSG00000095019 0 0 0 1
## FDR summary.logFC logFC.FALSE
## <numeric> <numeric> <numeric>
## ENSMUSG00000040204 3.72798e-46 3.65259 3.65259
## ENSMUSG00000017716 7.27422e-43 3.89331 3.89331
## ENSMUSG00000031004 1.46834e-37 3.30154 3.30154
## ENSMUSG00000069272 8.61905e-37 4.07536 4.07536
## ENSMUSG00000094777 3.86691e-35 3.55771 3.55771
## ... ... ... ...
## ENSMUSG00000095763 1 0 0
## ENSMUSG00000095523 1 0 0
## ENSMUSG00000095475 1 0 0
## ENSMUSG00000094855 1 0 0
## ENSMUSG00000095019 1 0 0chosen <- lapply(markers.up, function(x) {
x %>%
as.data.frame() %>%
tibble::rownames_to_column() %>%
filter(!grepl("^(Rps|Rpl)", Symbol)) %>%
filter(self.detected > 0.1) %>%
slice(n=1:10) %>%
dplyr::select(rowname)
})
chosen <- data.table::rbindlist(chosen, idcol="cluster")
chosen$cluster <- as.factor(chosen$cluster)
chosen$symbol <- rowData(sce)[chosen$rowname, "Symbol"]
m <- as.matrix(logcounts(sce)[chosen$rowname,])
n_cluster <- length(unique(colLabels(sce)))
colors = as.character(paletteer::paletteer_c("viridis::cividis", n=n_cluster))
names(colors) = unique(colLabels(sce))
col_fun = circlize::colorRamp2(seq(-2,2, length.out = 9),
rev(RColorBrewer::brewer.pal(name = "RdBu", n = 9)))
Heatmap(t(scale(t(m))),
use_raster = TRUE,
cluster_rows = FALSE,
cluster_columns = FALSE,
column_split = colLabels(sce),
row_split = chosen$cluster,
row_labels = chosen$symbol,
row_names_gp = grid::gpar(fontsize = 8),
top_annotation = HeatmapAnnotation(
cluster=colLabels(sce),
col = list(
cluster = cluster_colors
)
),
name = "z-score",
col = col_fun,
show_column_names = FALSE)## 'magick' package is suggested to install to give better rasterization.
##
## Set `ht_opt$message = FALSE` to turn off this message.
Annotate with Hudson
# create a bulk reference dds object
library(tximeta)
library(BiocParallel)
library(SingleR)
anno <- readr::read_csv(file = "reference_hudson/20200410_Hudson_Anno.csv")## Registered S3 method overwritten by 'cli':
## method from
## print.boxx spatstat.geom
## Rows: 24 Columns: 7## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (6): SRR_no, Subset_time, Subset, Dataset, PMID, Cells
## dbl (1): day##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.anno$files <- file.path("reference_hudson/RNA/", anno$SRR_no, "quant.sf")
anno$names <- anno$SRR_no
se <- tximeta(anno)## importing quantifications## reading in files with read_tsv## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## found matching transcriptome:
## [ GENCODE - Mus musculus - release M24 ]
## loading existing TxDb created: 2023-08-17 05:33:26
## Loading required package: GenomicFeatures
## Loading required package: AnnotationDbi
##
## Attaching package: 'AnnotationDbi'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## loading existing transcript ranges created: 2023-08-17 05:33:27
## fetching genome info for GENCODEgse <- summarizeToGene(se)## loading existing TxDb created: 2023-08-17 05:33:26
## obtaining transcript-to-gene mapping from database
## loading existing gene ranges created: 2023-08-17 05:33:31
## summarizing abundance
## summarizing counts
## summarizing lengthrownames(gse) = gsub("\\..*$", "", rownames(gse))
pred <- SingleR(test = sce, ref = gse,
labels = gse$Subset,
assay.type.test=1,
assay.type.ref=1,
BPPARAM = MulticoreParam())
plotScoreHeatmap(pred,
annotation_col=as.data.frame(colData(sce)[,"label",drop=FALSE]))
sce$hudson_labels <- pred$pruned.labelscd8 <- sce$label %in% c(1,3,4,5,8)
tableau10medium = c("#729ECE", "#FF9E4A", "#67BF5C", "#ED665D",
"#AD8BC9", "#A8786E", "#ED97CA", "#A2A2A2",
"#CDCC5D", "#6DCCDA")
p_husdon_1 <- plotUMAP(sce[,cd8], colour_by = "hudson_labels") +
scale_color_manual(values = c(
`Transitory` = tableau10medium[1],
`Stem-like` = tableau10medium[3],
`Exhausted` = tableau10medium[4],
`Naive` = tableau10medium[8]
))## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.p_husdon_1
plotUMAP(sce[,cd8], colour_by = "hudson_labels",
point_alpha=0.3) +
scale_color_manual(values = c(
`Transitory` = tableau10medium[1],
`Stem-like` = tableau10medium[3],
`Exhausted` = tableau10medium[4],
`Naive` = tableau10medium[8]
)) +
geom_density_2d(aes(colour=colour_by),
contour_var = "ndensity")## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plotTableBarplot(x = sce[,cd8]$Sample,
y = sce[,cd8]$hudson_labels,
label_x = "Sample",
label_y = "Hudson label") +
scale_fill_manual(values = c(
`Transitory` = tableau10medium[1],
`Stem-like` = tableau10medium[3],
`Exhausted` = tableau10medium[4],
`Naive` = tableau10medium[8]
))## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
T cell function
Pathway analysis with AUCell
AUCell
library(AUCell)
cell_rankings <- AUCell::AUCell_buildRankings(logcounts(sce), nCores = 1)## Quantiles for the number of genes detected by cell:
## (Non-detected genes are shuffled at the end of the ranking. Keep it in mind when choosing the threshold for calculating the AUC).
## min 1% 5% 10% 50% 100%
## 339 369 408 449 672 4163TRM signature
Milner = read.csv(file = "gene_lists/2021-04-22 Core Trm signature_Milner et al Nature 2017_vIL_corrected.txt",
header = FALSE)[[1]]
trm <- list(
Milner = rownames(sce)[match(Milner, rowData(sce)$Symbol)]
)
cell_AUC_trm <- AUCell_calcAUC(trm, cell_rankings, aucMaxRank=nrow(cell_rankings)*0.05,
verbose=FALSE, nCores=10)
## add to colData
colData(sce) <- cbind(colData(sce), t(getAUC(cell_AUC_trm)))
plotColData(sce, x="Sample", y="Milner", colour_by = "Sample") +
geom_boxplot(outlier.shape = NA,
width = .2,
alpha= .4) +
ggh4x::force_panelsizes(rows = unit(4, "in"),
cols = unit(2, "in"))
Exhaustion signature
exh <- read.csv("gene_lists/exhaustion_gene_lists.csv")
exh <- list(
Wherry = na.omit(rownames(sce)[match(na.omit(exh$Wherry), rowData(sce)$Symbol)])
)
cell_AUC_exh <- AUCell_calcAUC(exh, cell_rankings, aucMaxRank=nrow(cell_rankings)*0.05,
verbose=FALSE, nCores=10)
## add to colData
colData(sce) <- cbind(colData(sce), t(getAUC(cell_AUC_exh)))
plotColData(sce, x="Sample", y="Wherry", colour_by = "Sample") +
geom_boxplot(outlier.shape = NA,
width = .2,
alpha= .4) +
ggh4x::force_panelsizes(rows = unit(4, "in"),
cols = unit(2, "in"))
Synergistic effects of ICB and Propranolol
CD8 T cells
To test that, we compare ICB vs IgG2a, Propranolol vs IgG2a and Propranolol_ICB vs IgG2a.
markers <- findMarkers(sce[,cd8], direction="up", groups=sce[,cd8]$Sample, full.stats=TRUE,
test.type="wilcox")
markers_up <- list(
ICB = rownames(subset(markers$ICB$stats.IgG2A, log.FDR < log(0.05))),
Propranolol = rownames(subset(markers$Propranolol$stats.IgG2A, log.FDR < log(0.05))),
Propranolol_ICB = rownames(subset(markers$Propranolol_ICB$stats.IgG2A, log.FDR < log(0.05)))
)
list_to_upset_dataframe <- function(l) {
u = unique(unlist(l))
df = data.frame(gene = u)
for(n in names(l))
df[[n]] = u %in% l[[n]]
df
}
venn_df = list_to_upset_dataframe(markers_up)
ensembl2gene = as.data.frame(rowData(sce))[,1:2]
venn_df %>%
dplyr::left_join(ensembl2gene, by = c("gene" = "ID")) %>%
arrange(desc(Propranolol_ICB), desc(Propranolol), desc(ICB)) %>%
dplyr::select(Symbol, everything()) %>%
DT::datatable(
rownames = FALSE,
filter = "top",
extensions = 'Buttons',
options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel')
)
)Session Info
sessionInfo()## R version 4.1.3 (2022-03-10)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Fedora Linux 36 (Container Image)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libflexiblas.so.3.3
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices datasets
## [8] utils methods base
##
## other attached packages:
## [1] AUCell_1.14.0 GenomicFeatures_1.44.2
## [3] AnnotationDbi_1.54.1 SingleR_1.6.1
## [5] tximeta_1.10.0 dplyr_1.0.7
## [7] ComplexHeatmap_2.8.0 patchwork_1.1.1
## [9] magieR_0.1.0 BiocParallel_1.26.2
## [11] scran_1.20.1 scater_1.20.1
## [13] ggplot2_3.3.5 scuttle_1.2.1
## [15] SingleCellExperiment_1.14.1 SummarizedExperiment_1.22.0
## [17] Biobase_2.52.0 GenomicRanges_1.44.0
## [19] GenomeInfoDb_1.28.2 IRanges_2.26.0
## [21] S4Vectors_0.30.0 BiocGenerics_0.38.0
## [23] MatrixGenerics_1.4.3 matrixStats_0.60.1
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 rtracklayer_1.52.1
## [3] scattermore_0.7 ggthemes_4.2.4
## [5] R.methodsS3_1.8.1 SeuratObject_4.0.2
## [7] tidyr_1.1.3 bit64_4.0.5
## [9] knitr_1.33 R.utils_2.10.1
## [11] irlba_2.3.3 DelayedArray_0.18.0
## [13] data.table_1.14.0 rpart_4.1-15
## [15] AnnotationFilter_1.16.0 KEGGREST_1.32.0
## [17] RCurl_1.98-1.4 doParallel_1.0.16
## [19] generics_0.1.0 ScaledMatrix_1.0.0
## [21] cowplot_1.1.1 RSQLite_2.2.8
## [23] RANN_2.6.1 future_1.22.1
## [25] tzdb_0.1.2 bit_4.0.4
## [27] xml2_1.3.2 spatstat.data_2.1-0
## [29] httpuv_1.6.2 isoband_0.2.5
## [31] assertthat_0.2.1 viridis_0.6.1
## [33] tximport_1.20.0 xfun_0.25
## [35] hms_1.1.0 jquerylib_0.1.4
## [37] evaluate_0.14 promises_1.2.0.1
## [39] restfulr_0.0.13 fansi_0.5.0
## [41] progress_1.2.2 dbplyr_2.1.1
## [43] igraph_1.2.6.9118 DBI_1.1.1
## [45] htmlwidgets_1.5.3 spatstat.geom_2.2-2
## [47] paletteer_1.4.0 purrr_0.3.4
## [49] ellipsis_0.3.2 crosstalk_1.1.1
## [51] RSpectra_0.16-0 bookdown_0.23
## [53] annotate_1.70.0 prismatic_1.0.0
## [55] biomaRt_2.48.3 deldir_0.2-10
## [57] sparseMatrixStats_1.4.2 vctrs_0.3.8
## [59] ensembldb_2.16.4 here_1.0.1
## [61] Cairo_1.5-12.2 ROCR_1.0-11
## [63] abind_1.4-5 cachem_1.0.6
## [65] withr_2.4.2 ggh4x_0.2.1
## [67] vroom_1.5.4 sctransform_0.3.2
## [69] GenomicAlignments_1.28.0 prettyunits_1.1.1
## [71] goftest_1.2-2 cluster_2.1.2
## [73] dir.expiry_1.0.0 lazyeval_0.2.2
## [75] crayon_1.4.1 basilisk.utils_1.4.0
## [77] edgeR_3.34.0 pkgconfig_2.0.3
## [79] labeling_0.4.2 ProtGenerics_1.24.0
## [81] nlme_3.1-152 vipor_0.4.5
## [83] rlang_0.4.11 globals_0.14.0
## [85] lifecycle_1.0.0 miniUI_0.1.1.1
## [87] filelock_1.0.2 BiocFileCache_2.0.0
## [89] rsvd_1.0.5 AnnotationHub_3.0.1
## [91] rprojroot_2.0.2 polyclip_1.10-0
## [93] lmtest_0.9-38 graph_1.70.0
## [95] Matrix_1.3-4 zoo_1.8-9
## [97] beeswarm_0.4.0 pheatmap_1.0.12
## [99] ggridges_0.5.3 GlobalOptions_0.1.2
## [101] png_0.1-7 viridisLite_0.4.0
## [103] rjson_0.2.20 bitops_1.0-7
## [105] R.oo_1.24.0 KernSmooth_2.23-20
## [107] Biostrings_2.60.2 blob_1.2.2
## [109] DelayedMatrixStats_1.14.3 shape_1.4.6
## [111] stringr_1.4.0 parallelly_1.27.0
## [113] readr_2.0.1 beachmat_2.8.1
## [115] scales_1.1.1 GSEABase_1.54.0
## [117] memoise_2.0.0 magrittr_2.0.1
## [119] plyr_1.8.6 ica_1.0-2
## [121] zlibbioc_1.38.0 compiler_4.1.3
## [123] BiocIO_1.2.0 dqrng_0.3.0
## [125] RColorBrewer_1.1-2 clue_0.3-59
## [127] fitdistrplus_1.1-5 cli_3.0.1
## [129] Rsamtools_2.8.0 XVector_0.32.0
## [131] listenv_0.8.0 pbapply_1.4-3
## [133] MASS_7.3-54 mgcv_1.8-36
## [135] tidyselect_1.1.1 stringi_1.7.4
## [137] highr_0.9 yaml_2.2.1
## [139] BiocSingular_1.8.1 locfit_1.5-9.4
## [141] ggrepel_0.9.1 sass_0.4.0
## [143] tools_4.1.3 future.apply_1.8.1
## [145] rstudioapi_0.13 circlize_0.4.13
## [147] bluster_1.2.1 foreach_1.5.1
## [149] metapod_1.0.0 gridExtra_2.3
## [151] rmdformats_1.0.2 farver_2.1.0
## [153] Rtsne_0.15 digest_0.6.27
## [155] BiocManager_1.30.16 shiny_1.6.0
## [157] Rcpp_1.0.7 BiocVersion_3.13.1
## [159] later_1.3.0 RcppAnnoy_0.0.19
## [161] httr_1.4.2 colorspace_2.0-2
## [163] XML_3.99-0.7 tensor_1.5
## [165] reticulate_1.24 splines_4.1.3
## [167] uwot_0.1.10 statmod_1.4.36
## [169] rematch2_2.1.2 spatstat.utils_2.2-0
## [171] basilisk_1.4.0 renv_0.15.4
## [173] plotly_4.9.4.1 xtable_1.8-4
## [175] jsonlite_1.7.2 R6_2.5.1
## [177] pillar_1.6.2 htmltools_0.5.2
## [179] mime_0.11 DT_0.18
## [181] glue_1.4.2 fastmap_1.1.0
## [183] BiocNeighbors_1.10.0 interactiveDisplayBase_1.30.0
## [185] codetools_0.2-18 utf8_1.2.2
## [187] lattice_0.20-44 bslib_0.2.5.1
## [189] spatstat.sparse_2.0-0 tibble_3.1.4
## [191] curl_4.3.2 ggbeeswarm_0.6.0
## [193] leiden_0.3.9 survival_3.2-11
## [195] limma_3.48.3 rmarkdown_2.10
## [197] munsell_0.5.0 GetoptLong_1.0.5
## [199] GenomeInfoDbData_1.2.6 iterators_1.0.13
## [201] reshape2_1.4.4 gtable_0.3.0
## [203] spatstat.core_2.3-0 Seurat_4.0.4